#### "Economic Inequality. Income and political participation in authoritarian regimes" ####
# authors: "Pelke, Lars"
# date: 2020-06-25
# written under "R version 3.6.0 (2020-06-22)"
# Robustness test vote models 

# Structure #
## 1. Income Deciles 
## 2. All the Ginis measure
## 3. Test: Compulsory voting
## 5. Multiple Imputation of Missings

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# set working directory

# loading packages

library(countrycode)
library(tidyverse)
library(viridis)
library(haven)
library(lme4)
library(sjstats)
library(texreg)
library(interplot)
library(broom.mixed)

###############################################################################################
###############################################################################################
#### 1. Income Deciles 
#### Supplementary Appendix D1 and D2

#### Import Data ####
wvsdata <- readRDS("data/wvsdata.rds") 

#### Reducing to authoritarian regimes ####

wvsdata <- wvsdata %>%
  filter(v2x_regime_amb<=4)

#### Analyze Regime with elections multiparty:
# at least "constrained. at least one real opposition party is 
# allowed to context but competition is highly constrained -legally or informally. 

wvsdata <- wvsdata %>%
  filter(v2elmulpar_ord >=2)

wvsdata <- wvsdata %>%
  mutate(country_year = stringr::str_c(country_name, year, sep=" "))

table(wvsdata$country_year, is.na(wvsdata$gini_mkt))
table(wvsdata$country_year, is.na(wvsdata$v2elvotbuy))

number_observations <- wvsdata %>%
  count(country_year)


#### Building variables ####

#Replacing missings in Dependent Variables

wvsdata <- wvsdata %>%
  filter(E264>=0 | E257>=0) %>%
  mutate(vote = ifelse(E264==1, 1,
                       ifelse(E264==2, 1,
                              ifelse(E257==1, 1, 0))))
table(wvsdata$vote)

#Repalcing missings in Independent Variables

wvsdata <- wvsdata %>%
  filter(X001>=0)%>%
  mutate(sex= X001)

wvsdata <- wvsdata %>%
  filter(X003>=0)%>%
  mutate(age= X003)

wvsdata <- wvsdata %>%
  filter(X025>=0)%>%
  mutate(education= X025)

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(employement_status= X028)

wvsdata <- wvsdata %>%
  filter(X047>=0)%>%
  mutate(income_deciles= X047)

wvsdata <- wvsdata %>%
  mutate(income_quintiles = ifelse(income_deciles==1, 1,
                                   ifelse(income_deciles==2, 1,
                                          ifelse(income_deciles==3, 2,
                                                 ifelse(income_deciles==4, 2,
                                                        ifelse(income_deciles==5, 3,
                                                               ifelse(income_deciles==6, 3,
                                                                      ifelse(income_deciles==7 , 4,
                                                                             ifelse(income_deciles==8, 4, 5)))))))))


wvsdata <- wvsdata %>%
  filter(X011>=0)%>%
  mutate(children= ifelse(X011>0, 1, 0))

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(unemployed= ifelse(X028==7, 1, 0))

table(wvsdata$vote)
# 11236 non vote, 31942 vote

number_observations2 <- wvsdata %>%
  count(country_year)

## Inspect which countries are delete because missing obserations at macro-level variables ##

setdiff(number_observations$country_year, number_observations2$country_year)
#  [1] "Albania 1998"         "Algeria 2002"         "Armenia 1997"         "Azerbaijan 1997"      "Bangladesh 2002"     
# [6] "Belarus 1990"         "Croatia 1996"         "Egypt 2001"           "Georgia 1996"         "Hong Kong 2005"      
# [11] "Hong Kong 2014"       "Iran 2000"            "Iran 2007"            "Jordan 2001"          "Jordan 2007"         
# [16] "Kyrgyzstan 2003"      "Mexico 1990"          "Morocco 2001"         "Nigeria 1990"         "Nigeria 1995"        
# [21] "Nigeria 2000"         "North Macedonia 1998" "North Macedonia 2001" "Pakistan 1997"        "Pakistan 2001"       
# [26] "Peru 1996"            "Poland 1989"          "Russia 1990"          "Russia 1995"          "Rwanda 2007"         
# [31] "Singapore 2002"       "South Africa 1982"    "South Africa 1990"    "South Korea 1982"     "Taiwan 1994"         
# [36] "Tanzania 2001"        "Uganda 2001"          "Zanzibar 2001"        "Zimbabwe 2001"      

###Data Manipulation of GINI measures, fill some missing value with last national observations

missing <- wvsdata %>% 
  group_by(country_year) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
missing
rm(missing)


wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Algeria 2014", 36.7, gini_mkt)) %>% # Algeria 2011
  mutate(gini_disp=ifelse(country_year=="Algeria 2014", 35, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Azerbaijan 2011", 38.3, gini_mkt)) %>% # Azerbaijan 2011
  mutate(gini_disp=ifelse(country_year=="Azerbaijan 2011", 30.3, gini_disp)) %>%
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 38.5, gini_mkt)) %>% # Lebanon 2012
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 36.5, gini_disp))

wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Haiti 2016", 57.6, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Haiti 2016", 53.7, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Nigeria 2012", 46.8, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Nigeria 2012", 44.1, gini_disp)) 


wvsdata <- wvsdata %>%
  drop_na(gini_mkt)

#### Previos democratic experience ####

table(wvsdata$democratic_expericene) 
table(wvsdata$democratic_expericene_all)

democratic_exp <- wvsdata %>%
  filter(democratic_expericene_bi ==1) %>%
  group_by(country_name) %>%
  count()
# authoritarian regime that has previous democratic experience: Armenia (1990-1994), Belarus (1991-1996), Palestine (2004-2006), Russia (1992, 1993), 
# Thailand (1998-2005, 2012)

rm(democratic_exp)

#### Social Group Base Regime ####

table(wvsdata$regime_support_group_max) 

#### Country and Country-Year overview ####

country_year_observations <- wvsdata %>% 
  group_by(country_year) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
country_year_observations
rm(country_year_observations)

country_observations <- wvsdata %>% 
  group_by(country_name) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
country_observations
rm(country_observations)

count <- wvsdata %>%
  group_by(S025, income_quintiles) %>%
  tally() 
rm(count)


number_observations3 <- wvsdata %>%
  count(country_year)
setdiff(number_observations2$country_year, number_observations3$country_year)
# [1] "Kuwait 2014"  "Lebanon 2013" "Libya 2014"  

rm(number_observations3, number_observations2, number_observations)

##################################
#Group and grand mean centering
##################################

# Author: Zoltan Fazekas
centFUN <- function(x) {
  x - mean(x, na.rm = TRUE)
}

# Median Cenetering for factor variables
centFUN_median <- function(x) {
  x - median(x, na.rm = TRUE)
}

#### Data manipulation ####

wvsdata <- wvsdata %>%
  mutate(democratic_expericene_all = ifelse(is.na(democratic_expericene_all), 0, democratic_expericene_all), 
         democratic_expericene = ifelse(is.na(democratic_expericene), 0, democratic_expericene))


# Use the function to grand-mean center the macro-level predictor
wvsdata$gini_mktCGM <- centFUN(wvsdata$gini_mkt)
wvsdata$gini_dispCGM <- centFUN(wvsdata$gini_disp)
wvsdata$v2xel_frefairCGM <- centFUN(wvsdata$v2xel_frefair)
wvsdata$e_migdppclnCGM <- centFUN(wvsdata$e_migdppcln)
wvsdata$v2elvotbuyCGM <- centFUN(wvsdata$v2elvotbuy)
wvsdata$regime_support_group_maxCGM <- centFUN(wvsdata$regime_support_group_max)
wvsdata$democratic_expericene_allCGM <- centFUN(wvsdata$democratic_expericene_all)
wvsdata$democratic_expericeneCGM <- centFUN(wvsdata$democratic_expericene)


summary(wvsdata$e_migdppclnCGM)
summary(wvsdata$v2xel_frefairCGM)
summary(wvsdata$gini_mktCGM)
summary(wvsdata$gini_dispCGM)
summary(wvsdata$v2elvotbuyCGM)
summary(wvsdata$regime_support_group_maxCGM)
summary(wvsdata$democratic_expericeneCGM)
summary(wvsdata$democratic_expericene_allCGM)

#Use the fuction to group-mean center the individual predictors
wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(ageCWC = centFUN(age))
summary(wvsdata$ageCWC)

wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(educationCWC = centFUN(as.numeric(education)),
         income_quintilesCWC = centFUN(as.numeric(income_quintiles)),
         income_decilesCWC = centFUN(as.numeric(income_deciles)))

summary(wvsdata$educationCWC)
summary(wvsdata$income_quintilesCWC)
summary(wvsdata$income_decilesCWC)

#Rescaling predicator variables
summary(wvsdata$educationCWC)
summary(wvsdata$sex)
summary(wvsdata$ageCWC)
summary(wvsdata$income_quintilesCWC)
summary(wvsdata$S003)
wvsdata$ageCWC <- wvsdata$ageCWC/10

##################################
#Analysis: three-level models: Market Inequality
##################################

v1 <- glmer(vote ~ 1 + (1 | S025) + (1 | S003) ,
            data = wvsdata, 
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v1)
icc(v1, adjusted = TRUE)
icc(v1)


v2 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v2)
anova(v1, v2)

#macrolevel gini
v3 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_mktCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v3)
anova(v2, v3)

#macrolevel gini plus additional macro controls
v4 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_mktCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v4)
anova(v2, v4)

# model 4 plus random slopes
v5 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_mktCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v5)
anova(v2, v5)

# model 5 plus interactions 
v6 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC +
              gini_mktCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_decilesCWC*gini_mktCGM + v2elvotbuyCGM*income_decilesCWC +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v6)
anova(v5, v6)


# model 6 but with interaction of freedom/fariness and income level and without gini
v7 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC +
              v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_decilesCWC*v2xel_frefairCGM +  v2elvotbuyCGM*income_decilesCWC +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v7)
anova(v6, v7)

##Interplot relationship between income and inequality

library(interplot)
plot_Interaction_vote_mkt <- interplot(v6, 
                                       var1 = "gini_mktCGM", 
                                       var2 = "income_decilesCWC", 
                                       ci = 0.95, 
                                       hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of inequality on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_mkt
ggsave("output/D1_interaction_inequality.pdf")


plot_Interaction_vote_buying <- interplot(v6, 
                                          var1 = "v2elvotbuyCGM", 
                                          var2 = "income_decilesCWC", 
                                          ci = 0.95, 
                                          hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of vote buying on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_buying
ggsave("output/D1_interaction_votebuying.pdf")


plot_Interaction_vote_cleanelections <- interplot(v7, 
                                                  var1 = "v2xel_frefairCGM", 
                                                  var2 = "income_decilesCWC", 
                                                  ci = 0.95, 
                                                  hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of clean elections on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_cleanelections
ggsave("output/D1_interaction_cleanelections.pdf")


######################################################
#Building Logg Odds ratios of Models

library(dotwhisker)
library(broom.mixed)


v1_df <- tidy(v1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v2_df <- tidy(v2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v3_df <- tidy(v3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v4_df <- tidy(v4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v5_df <- tidy(v5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v6_df <- tidy(v6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v7_df <- tidy(v7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
v1_df <- v1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v2_df <- v2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))
v3_df <- v3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v4_df <- v4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v5_df <- v5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v6_df <- v6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v7_df <- v7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

#Log Odds Plot Effects Plot

v_models <- rbind(v2_df, v3_df, v4_df, v5_df, v6_df, v7_df)

v_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors", "Market Inequality (centered)", "Regime Support Group (centered)"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * Fairness election"))

plot_main_vote <- {dwplot(v_models) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Voting") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.75, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())} %>% 
  add_brackets(v_brackets)

plot_main_vote
ggsave("output/D1_dot_whisker.pdf", width = 22 , height = 25 , units = c("cm"))


#### texreg regression tables main models ####

texreg(list(v1, v2, v3, v4, v5, v6, v7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income deciles (centered)",
                             "Inequality (centered)",
                             "Fairness elections (centered)",
                             "GDP pc log (centered)",
                             "Vote Buying (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * Fairness election",
                             "Income * Vote Buying"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Main Models Voting",
       fontsize = "scriptsize",
       label = "Table:D1", 
       longtable =  TRUE)

###############################
##Using Disposable Inequality##
###############################

d1 <- glmer(vote ~ 1 + (1 | S025) + (1 | S003) ,
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d1)

d2 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d2)
anova(d1, d2)

#macrolevel gini
d3 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_dispCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d3)
anova(d2, d3)

#macrolevel gini plus addtional macro controls
d4 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_dispCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d4)
anova(d2, d4)

# model 4 plus random slopes
d5 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_dispCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d5)
anova(d2, d5)

# model 5 plus random slopes
d6 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed+ income_decilesCWC +
              gini_dispCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_decilesCWC*gini_dispCGM +  income_decilesCWC*v2elvotbuyCGM +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d6)
anova(d5, d6)

# model 5 without Inequality and interaction
d7 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed+ income_decilesCWC +
              v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_decilesCWC*v2xel_frefairCGM +  income_decilesCWC*v2elvotbuyCGM +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d7)
anova(d5, d7)

##Interplot relationship between income and inequality

plot_Interaction_vote_disp <- interplot(d6, 
                                        var1 = "gini_dispCGM", 
                                        var2 = "income_decilesCWC", 
                                        ci = 0.95, 
                                        hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of disposable inequality on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_disp
ggsave("output/D2_interaction_inequality.pdf")


plot_Interaction_vote_disp_vote_buy <- interplot(d6, 
                                                 var1 = "v2elvotbuyCGM", 
                                                 var2 = "income_decilesCWC", 
                                                 ci = 0.95, 
                                                 hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of vote buying on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_disp_vote_buy
ggsave("output/D2_interaction_votebuying.pdf")


plot_Interaction_vote_disp_cleanelections <- interplot(d7, 
                                                       var1 = "v2xel_frefairCGM", 
                                                       var2 = "income_decilesCWC", 
                                                       ci = 0.95, 
                                                       hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of clean elections on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_disp_cleanelections
ggsave("output/D2_interaction_cleanelelections.pdf")


#### Plotting Dot-Whisker Plot main Models Vote ####

d1_df <- tidy(d1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d2_df <- tidy(d2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d3_df <- tidy(d3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d4_df <- tidy(d4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d5_df <- tidy(d5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d6_df <- tidy(d6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d7_df <- tidy(d7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
d1_df <- d1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d2_df <- d2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d3_df <- d3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d4_df <- d4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))


d5_df <- d5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))


d6_df <- d6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d7_df <- d7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))


#Log Odds Plot Effects Plot

v_models <- rbind(d2_df, d3_df, d4_df, d5_df, d6_df, d7_df)

v_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors", "Disp. Inequality (centered)", "Regime Support Group (centered)"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * Fairness election"))

plot_main_vote_disp <- {dwplot(v_models) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Voting") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.85, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())} %>% 
  add_brackets(v_brackets)

plot_main_vote_disp
ggsave("output/D2_dot_whisker.pdf", width = 22 , height = 25 , units = c("cm"))


#### Texreg regression table ####

texreg(list(d1, d2, d3, d4, d5, d6, d7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income deciles (centered)",
                             "Inequality (centered)",
                             "Fairness elections (centered)",
                             "GDP pc log (centered)",
                             "Vote Buying (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * Fairness election",
                             "Income * Vote Buying"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Main Models Voting (Disp. Inequality)",
       fontsize = "scriptsize",
       label = "Table:D2", 
       longtable = TRUE)

###############################################################################################
###############################################################################################
#### 2. All the Ginis Data 
#### Supplementary Appendix D3


# clear workspace
rm(list=ls())

# set working directory
setwd("M:/Promotion/Paper/WVS Paper Inequality and Voting in Autocracies/calculations_new")

# loading packages

library(countrycode)
library(tidyverse)
library(viridis)
library(haven)
library(lme4)
library(sjstats)
library(texreg)

#### Import Data ####
wvsdata <- readRDS("data/wvsdata.rds") 

#### Import all the Ginis measure by Vdem version 8.0 ####

vdem8 <- read_csv("data/vdem_8/V-Dem-CY+Others-v8.csv")
vdem8$cowcode <- vdem8$COWcode

vdem8 <- vdem8 %>%
  dplyr::select(country_id, year, e_peginiwi)

#### Merge All the ginis with WVSdata

wvsdata <- wvsdata %>%
  left_join(vdem8, c("country_id", "year"))

rm(vdem8)

#### Reducing to authoritarian regimes ####

wvsdata <- wvsdata %>%
  filter(v2x_regime_amb<=4)

#### Analyze Regime with elections multiparty:
# at least "constrained. at least one real opposition party is 
# allowed to context but competition is highly constrained -legally or informally. 

wvsdata <- wvsdata %>%
  filter(v2elmulpar_ord >=2)

wvsdata <- wvsdata %>%
  mutate(country_year = stringr::str_c(country_name, year, sep=" "))

table(wvsdata$country_year, is.na(wvsdata$e_peginiwi))
table(wvsdata$country_year, is.na(wvsdata$v2elvotbuy))

number_observations <- wvsdata %>%
  count(country_year)


#### Building variables ####

#Replacing missings in Dependent Variables

wvsdata <- wvsdata %>%
  filter(E264>=0 | E257>=0) %>%
  mutate(vote = ifelse(E264==1, 1,
                       ifelse(E264==2, 1,
                              ifelse(E257==1, 1, 0))))
table(wvsdata$vote)

#Repalcing missings in Independent Variables

wvsdata <- wvsdata %>%
  filter(X001>=0)%>%
  mutate(sex= X001)

wvsdata <- wvsdata %>%
  filter(X003>=0)%>%
  mutate(age= X003)

wvsdata <- wvsdata %>%
  filter(X025>=0)%>%
  mutate(education= X025)

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(employement_status= X028)

wvsdata <- wvsdata %>%
  filter(X047>=0)%>%
  mutate(income_deciles= X047)

wvsdata <- wvsdata %>%
  mutate(income_quintiles = ifelse(income_deciles==1, 1,
                                   ifelse(income_deciles==2, 1,
                                          ifelse(income_deciles==3, 2,
                                                 ifelse(income_deciles==4, 2,
                                                        ifelse(income_deciles==5, 3,
                                                               ifelse(income_deciles==6, 3,
                                                                      ifelse(income_deciles==7 , 4,
                                                                             ifelse(income_deciles==8, 4, 5)))))))))


wvsdata <- wvsdata %>%
  filter(X011>=0)%>%
  mutate(children= ifelse(X011>0, 1, 0))

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(unemployed= ifelse(X028==7, 1, 0))

table(wvsdata$vote)
# 11236 non vote, 31942 vote

###Data Manipulation of All the Ginis measures

missing <- wvsdata %>% 
  group_by(country_year) %>% 
  summarise(sumNA = sum(is.na(e_peginiwi)))
missing
rm(missing)

wvsdata <- wvsdata %>%
  drop_na(e_peginiwi)

#### Previos democratic experience ####

table(wvsdata$democratic_expericene) 
table(wvsdata$democratic_expericene_all)

democratic_exp <- wvsdata %>%
  filter(democratic_expericene_bi ==1) %>%
  group_by(country_name) %>%
  count()
# authoritarian regime that has previous democratic experience: Armenia (1990-1994), Belarus (1991-1996), Palestine (2004-2006), Russia (1992, 1993), 
# Thailand (1998-2005, 2012)

rm(democratic_exp)

#### Social Group Base Regime ####

table(wvsdata$regime_support_group_max) 


##################################
#Group and grand mean centering
##################################

# Author: Zoltan Fazekas
centFUN <- function(x) {
  x - mean(x, na.rm = TRUE)
}

# Median Cenetering for factor variables
centFUN_median <- function(x) {
  x - median(x, na.rm = TRUE)
}

#### Data manipulation ####

wvsdata <- wvsdata %>%
  mutate(democratic_expericene_all = ifelse(is.na(democratic_expericene_all), 0, democratic_expericene_all), 
         democratic_expericene = ifelse(is.na(democratic_expericene), 0, democratic_expericene))

# Use the function to grand-mean center the macro-level predictor
wvsdata$e_peginiwiCGM <- centFUN(wvsdata$e_peginiwi)
wvsdata$v2xel_frefairCGM <- centFUN(wvsdata$v2xel_frefair)
wvsdata$e_migdppclnCGM <- centFUN(wvsdata$e_migdppcln)
wvsdata$v2elvotbuyCGM <- centFUN(wvsdata$v2elvotbuy)
wvsdata$regime_support_group_maxCGM <- centFUN(wvsdata$regime_support_group_max)
wvsdata$democratic_expericene_allCGM <- centFUN(wvsdata$democratic_expericene_all)
wvsdata$democratic_expericeneCGM <- centFUN(wvsdata$democratic_expericene)

summary(wvsdata$e_migdppclnCGM)
summary(wvsdata$v2xel_frefairCGM)
summary(wvsdata$e_peginiwiCGM)
summary(wvsdata$v2elvotbuyCGM)
summary(wvsdata$regime_support_group_maxCGM)
summary(wvsdata$democratic_expericeneCGM)
summary(wvsdata$democratic_expericene_allCGM)

#Use the fuction to group-mean center the individual predictors
wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(ageCWC = centFUN(age))
summary(wvsdata$ageCWC)

wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(educationCWC = centFUN(as.numeric(education)),
         income_quintilesCWC = centFUN(as.numeric(income_quintiles)),
         income_decilesCWC = centFUN(as.numeric(income_deciles)))

summary(wvsdata$educationCWC)
summary(wvsdata$income_quintilesCWC)
summary(wvsdata$income_decilesCWC)

#Rescaling predicator variables
wvsdata$ageCWC <- wvsdata$ageCWC/10

###############################
##Using UNU-WIDER DATA BY VDEM 8##
###############################

d1 <- glmer(vote ~ 1 + (1 | S025) + (1 | S003) ,
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d1)

d2 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d2)
anova(d1, d2)

#macrolevel gini
d3 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              e_peginiwiCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d3)
anova(d2, d3)

#macrolevel gini plus addtional macro controls
d4 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              e_peginiwiCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d4)
anova(d2, d4)

# model 4 plus random slopes
d5 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              e_peginiwiCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d5)
anova(d2, d5)

# model 5 plus random slopes
d6 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed+ income_decilesCWC +
              e_peginiwiCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_decilesCWC*e_peginiwiCGM +  income_decilesCWC*v2elvotbuyCGM +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d6)
anova(d5, d6)

# model 5 without Inequality and interaction
d7 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed+ income_decilesCWC +
              v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_decilesCWC*v2xel_frefairCGM +  income_decilesCWC*v2elvotbuyCGM +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d7)
anova(d5, d7)

##Interplot relationship between income and inequality

plot_Interaction_vote_disp <- interplot(d6, 
                                        var1 = "e_peginiwiCGM", 
                                        var2 = "income_decilesCWC", 
                                        ci = 0.95, 
                                        hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of inequality on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_disp
ggsave("output/D3interactioninequality.pdf")


plot_Interaction_vote_disp_vote_buy <- interplot(d6, 
                                                 var1 = "v2elvotbuyCGM", 
                                                 var2 = "income_decilesCWC", 
                                                 ci = 0.95, 
                                                 hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of vote buying on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_disp_vote_buy
ggsave("output/D3interactionvotebuying.pdf")


plot_Interaction_vote_disp_cleanelections <- interplot(d7, 
                                                       var1 = "v2xel_frefairCGM", 
                                                       var2 = "income_decilesCWC", 
                                                       ci = 0.95, 
                                                       hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of clean elections on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_disp_cleanelections
ggsave("output/D3interactioncleanelelections.pdf")


#### Plotting Dot-Whisker Plot main Models Vote ####
library(dotwhisker)
library(broom.mixed)

d1_df <- tidy(d1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d2_df <- tidy(d2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d3_df <- tidy(d3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d4_df <- tidy(d4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d5_df <- tidy(d5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d6_df <- tidy(d6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d7_df <- tidy(d7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
d1_df <- d1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d2_df <- d2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))


d3_df <- d3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d4_df <- d4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))


d5_df <- d5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d6_df <- d6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d7_df <- d7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_decilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_decilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

#Log Odds Plot Effects Plot

v_models <- rbind(d2_df, d3_df, d4_df, d5_df, d6_df, d7_df)

v_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors", "Inequality (centered)", "Regime Support Group (centered)"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * Fairness election"))

plot_main_vote_disp <- {dwplot(v_models) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Voting") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.85, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())} %>% 
  add_brackets(v_brackets)

plot_main_vote_disp
ggsave("output/D3dotwhisker.pdf", width = 22 , height = 25 , units = c("cm"))


#### Texreg regression table ####

texreg(list(d1, d2, d3, d4, d5, d6, d7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income deciles (centered)",
                             "Inequality (centered)",
                             "Fairness elections (centered)",
                             "GDP pc log (centered)",
                             "Vote Buying (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * Fairness election",
                             "Income * Vote Buying"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Main Models Voting (Disp. Inequality)",
       fontsize = "scriptsize",
       label = "Table:D3", 
       longtable = TRUE)


###############################################################################################
###############################################################################################
#### 3. Test: Compulsory Voting
#### Supplementary Appendix D4


#### Import Data ####
wvsdata <- readRDS("data/wvsdata.rds") 


#### Reducing to authoritarian regimes ####

wvsdata <- wvsdata %>%
  filter(v2x_regime_amb<=4)

#### Analyze Regime with elections multiparty:
# at least "constrained. at least one real opposition party is 
# allowed to context but competition is highly constrained -legally or informally. 

wvsdata <- wvsdata %>%
  filter(v2elmulpar_ord >=2)

wvsdata <- wvsdata %>%
  mutate(country_year = stringr::str_c(country_name, year, sep=" "))

table(wvsdata$country_year, is.na(wvsdata$gini_mkt))
table(wvsdata$country_year, is.na(wvsdata$v2elvotbuy))

number_observations <- wvsdata %>%
  count(country_year)

#### Building variables ####

#Replacing missings in Dependent Variables

wvsdata <- wvsdata %>%
  filter(E264>=0 | E257>=0) %>%
  mutate(vote = ifelse(E264==1, 1,
                       ifelse(E264==2, 1,
                              ifelse(E257==1, 1, 0))))
table(wvsdata$vote)

#Repalcing missings in Independent Variables

wvsdata <- wvsdata %>%
  filter(X001>=0)%>%
  mutate(sex= X001)

wvsdata <- wvsdata %>%
  filter(X003>=0)%>%
  mutate(age= X003)

wvsdata <- wvsdata %>%
  filter(X025>=0)%>%
  mutate(education= X025)

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(employement_status= X028)

wvsdata <- wvsdata %>%
  filter(X047>=0)%>%
  mutate(income_deciles= X047)

wvsdata <- wvsdata %>%
  mutate(income_quintiles = ifelse(income_deciles==1, 1,
                                   ifelse(income_deciles==2, 1,
                                          ifelse(income_deciles==3, 2,
                                                 ifelse(income_deciles==4, 2,
                                                        ifelse(income_deciles==5, 3,
                                                               ifelse(income_deciles==6, 3,
                                                                      ifelse(income_deciles==7 , 4,
                                                                             ifelse(income_deciles==8, 4, 5)))))))))


wvsdata <- wvsdata %>%
  filter(X011>=0)%>%
  mutate(children= ifelse(X011>0, 1, 0))

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(unemployed= ifelse(X028==7, 1, 0))

#Constructing Complusory Voting Control Variable
wvsdata <- wvsdata %>%
  mutate(comp_vote = if_else(country_year=="Egypt 2001", 1,
                             if_else(country_year=="Egypt 2008", 1,
                             if_else(country_year=="Egypt 2012", 1,
                                     if_else(country_year=="Singapore 2002", 1,
                                             if_else(country_year=="Singapore 2012", 1,
                                             if_else(country_year=="Thailand 2007", 1,
                                                     if_else(country_year=="Thailand 2013", 1, 0))))))))

table(wvsdata$vote)
# 11236 non vote, 31942 vote

###Data Manipulation of GINI measures, fill some missing value with last national observations

missing <- wvsdata %>% 
  group_by(country_year) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
missing
rm(missing)


wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Algeria 2014", 36.7, gini_mkt)) %>% # Algeria 2011
  mutate(gini_disp=ifelse(country_year=="Algeria 2014", 35, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Azerbaijan 2011", 38.3, gini_mkt)) %>% # Azerbaijan 2011
  mutate(gini_disp=ifelse(country_year=="Azerbaijan 2011", 30.3, gini_disp)) %>%
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 38.5, gini_mkt)) %>% # Lebanon 2012
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 36.5, gini_disp))

wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Haiti 2016", 57.6, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Haiti 2016", 53.7, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Nigeria 2012", 46.8, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Nigeria 2012", 44.1, gini_disp)) 


wvsdata <- wvsdata %>%
  drop_na(gini_mkt)

#### Previos democratic experience ####

table(wvsdata$democratic_expericene) 
table(wvsdata$democratic_expericene_all)

democratic_exp <- wvsdata %>%
  filter(democratic_expericene_bi ==1) %>%
  group_by(country_name) %>%
  count()
# authoritarian regime that has previous democratic experience: Armenia (1990-1994), Belarus (1991-1996), Palestine (2004-2006), Russia (1992, 1993), 
# Thailand (1998-2005, 2012)

rm(democratic_exp)

#### Social Group Base Regime ####

table(wvsdata$regime_support_group_max) 


##################################
#Group and grand mean centering
##################################

# Author: Zoltan Fazekas
centFUN <- function(x) {
  x - mean(x, na.rm = TRUE)
}

# Median Cenetering for factor variables
centFUN_median <- function(x) {
  x - median(x, na.rm = TRUE)
}

#### Data manipulation ####

wvsdata <- wvsdata %>%
  mutate(democratic_expericene_all = ifelse(is.na(democratic_expericene_all), 0, democratic_expericene_all), 
         democratic_expericene = ifelse(is.na(democratic_expericene), 0, democratic_expericene))


# Use the function to grand-mean center the macro-level predictor
wvsdata$gini_mktCGM <- centFUN(wvsdata$gini_mkt)
wvsdata$gini_dispCGM <- centFUN(wvsdata$gini_disp)
wvsdata$v2xel_frefairCGM <- centFUN(wvsdata$v2xel_frefair)
wvsdata$e_migdppclnCGM <- centFUN(wvsdata$e_migdppcln)
wvsdata$v2elvotbuyCGM <- centFUN(wvsdata$v2elvotbuy)
wvsdata$regime_support_group_maxCGM <- centFUN(wvsdata$regime_support_group_max)
wvsdata$democratic_expericene_allCGM <- centFUN(wvsdata$democratic_expericene_all)
wvsdata$democratic_expericeneCGM <- centFUN(wvsdata$democratic_expericene)


summary(wvsdata$e_migdppclnCGM)
summary(wvsdata$v2xel_frefairCGM)
summary(wvsdata$gini_mktCGM)
summary(wvsdata$gini_dispCGM)
summary(wvsdata$v2elvotbuyCGM)
summary(wvsdata$regime_support_group_maxCGM)
summary(wvsdata$democratic_expericeneCGM)
summary(wvsdata$democratic_expericene_allCGM)


#Use the fuction to group-mean center the individual predictors
wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(ageCWC = centFUN(age))
summary(wvsdata$ageCWC)

wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(educationCWC = centFUN(as.numeric(education)),
         income_quintilesCWC = centFUN(as.numeric(income_quintiles)),
         income_decilesCWC = centFUN(as.numeric(income_deciles)))

summary(wvsdata$educationCWC)
summary(wvsdata$income_quintilesCWC)
summary(wvsdata$income_decilesCWC)

#Rescaling predicator variables
summary(wvsdata$educationCWC)
summary(wvsdata$sex)
summary(wvsdata$ageCWC)
summary(wvsdata$income_quintilesCWC)
summary(wvsdata$S003)
wvsdata$ageCWC <- wvsdata$ageCWC/10

##################################
#Analysis: three-level models: Market Inequality
##################################

#macrolevel gini plus addtional macro controls
v4 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_mktCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + comp_vote +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v4)
anova(v2, v4)

# model 4 plus random slopes
v5 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_mktCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + comp_vote +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v5)
anova(v2, v5)

# model 5 plus interactions 
v6 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              gini_mktCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + comp_vote +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*gini_mktCGM + v2elvotbuyCGM*income_quintilesCWC +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v6)
anova(v5, v6)


# model 6 but with interaction of freedom/fariness and income level and without gini
v7 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + comp_vote +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*v2xel_frefairCGM +  v2elvotbuyCGM*income_quintilesCWC +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v7)
anova(v6, v7)

##Interplot relationship between income and inequality

library(interplot)
plot_Interaction_vote_mkt <- interplot(v6, 
                                       var1 = "gini_mktCGM", 
                                       var2 = "income_quintilesCWC", 
                                       ci = 0.95, 
                                       hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of inequality on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_mkt
ggsave("output/D5interactioninequality.pdf")


plot_Interaction_vote_buying <- interplot(v6, 
                                          var1 = "v2elvotbuyCGM", 
                                          var2 = "income_quintilesCWC", 
                                          ci = 0.95, 
                                          hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of vote buying on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_buying
ggsave("output/D5interactionvotebuying.pdf")

plot_Interaction_vote_cleanelections <- interplot(v7, 
                                                  var1 = "v2xel_frefairCGM", 
                                                  var2 = "income_quintilesCWC", 
                                                  ci = 0.95, 
                                                  hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of clean elections on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_cleanelections
ggsave("output/D5interactioncleanelections.pdf")

######################################################
#Building Logg Odds ratios of Models

library(dotwhisker)
library(broom.mixed)

v4_df <- tidy(v4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v5_df <- tidy(v5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v6_df <- tidy(v6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v7_df <- tidy(v7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations

v4_df <- v4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       comp_vote = "Compulsory Voting",
                       democratic_expericeneCGM = "Previos Dem. Experience", 
                       regime_support_group_maxCGM = "Regime Support Group",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v5_df <- v5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       comp_vote = "Compulsory Voting",
                       democratic_expericeneCGM = "Previos Dem. Experience", 
                       regime_support_group_maxCGM = "Regime Support Group",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v6_df <- v6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       comp_vote = "Compulsory Voting",
                       democratic_expericeneCGM = "Previos Dem. Experience", 
                       regime_support_group_maxCGM = "Regime Support Group",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v7_df <- v7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       comp_vote = "Compulsory Voting",
                       democratic_expericeneCGM = "Previos Dem. Experience", 
                       regime_support_group_maxCGM = "Regime Support Group",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

#Log Odds Plot Effects Plot

v_models <- rbind(v4_df, v5_df, v6_df, v7_df)

v_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors", "Market Inequality (centered)", "Regime Support Group"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * Fairness election"))

plot_main_vote <- {dwplot(v_models) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Voting") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.85, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())} %>% 
  add_brackets(v_brackets)

plot_main_vote
ggsave("output/D5dotwhisker.pdf", width = 22 , height = 25 , units = c("cm"))


#### texreg regression tables main models ####

texreg(list(v4, v5, v6, v7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income quintiles (centered)",
                             "Inequality (centered)",
                             "Fairness elections (centered)",
                             "GDP pc log (centered)",
                             "Vote Buying (centered)",
                             "Complusory Voting",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * Fairness election",
                             "Income * Vote Buying"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Main Models Voting",
       fontsize = "scriptsize",
       label = "Table:D4", 
       longtable = TRUE)
